Work probability distribution in systems driven out of equilibrium 
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We derive the difFerential equation describing the time evolution of the work probability distri- 
bution function of a stochastic system which is driven out of equilibrium by the manipulation of a 
parameter. We consider both systems described by their microscopic state or by a collective variable 
which identifies a quasiequilibrium state. We show that the work probability distribution can be 
represented by a path integral, which is dominated by "classical" paths in the large system size 
limit. We compare these results with simulated manipulation of mean-field systems. We discuss the 
range of applicability of the Jarzynski equality for evaluating the system free energy using these out- 
of-equilibrium manipulations. Large fluctuations in the work and the shape of the work distribution 
tails are also discussed. 
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Recent improvements in micromanipulation techniques 
have made it possible to observe experimentaUy work 
fluctuations and to measure the probability distribution 
of the work exerted on a system subject to external ma- 
nipulation. In particular, the probability distribution 
of the work has been measured in RNA pulling exper- 
iments ^ and for micrometer-sized colloidal particles 
dragged through a fluid |^ . Usually, because of technical 
limitations, this class of experiments is characterized by 
time scales much faster than the typical system relax- 
ation time. This hinders the possibility to perform the 
experiments in quasistatic conditions and thus to obtain 
direct measurements of the system thermodynamic state 
variables. The importance of the knowledge of work dis- 
tributions in such experiments resides in the fact that 
one can evaluate the free energy difference between the 
final and the initial state of the system by exploiting the 
Jarzynski equality (JE) ||, ||| 

(e-^^) = e-'^^^. (1) 

According to previous works ||^, a precise knowledge of 
the tails in the distributions provides information on how 
many experiments are needed in order to evaluate cor- 
rectly the free energy difference of a system using non- 
equilibrium experiments. Thus, a priori estimates of 
P{W) are in principle needed, to evaluate the actual use- 
fulness of this approach. 

In two recent works || , we introduced and discussed 
a differential equation describing the time evolution of 
the probability distribution of the work done on a system 
by manipulating an external field (force) /i, according to 
a given protocol /i(t). In particular, in ref. |Q, we con- 
sidered the case of a system characterized by a discrete 
phase space, while in ref. |^ we considered a mean field 
system characterized by a generic equilibrium free energy 
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The aim of this paper is to extend those works, by 
exploiting an approach due to Felix Ritort |^ . In partic- 
ular, we first derive explicitly the differential equations 
governing the time evolution of P{W,t). We then de- 
rive an expression of the work probability distribution of 
a system described by a collective variable, on the hy- 
pothesis that, during the manipulation, the system finds 
itself in a quasiequilibrium state constrained by the value 
of that coordinate. We solve the resulting equation by 
path integrals and show that, in the limit of large system 
size, the path integral is dominated by the classical path 
which satisfy canonical equations of motion, and suitable 
boundary conditions. The expression for the probability 
distribution function follows straightforwardly. We high- 
light the analogy between the path functionals obtained 
in this way and classical thermodynamics. We apply the 
obtained results to some simple systems, and we explore 
in particular the possibility of the existence of exponen- 
tial tails in the work probability distribution: such tails 
are related, via the thermodynamic analogy, to phase 
transitions in the path distribution. We show that, con- 
trary to what was conjectured in ref. |^ on the basis of 
numerical evidence, such tails are not present in a para- 
magnet, or in a ferromagnet above the critical tempera- 
ture, but are present in a mean-field ferromagnet below 
the critical temperature, provided the manipulating pro- 
tocol is fast enough. The implications of our results are 
further discussed. 



I. PROBABILITY DISTRIBUTION OF THE 
WORK FOR THE MICROSCOPIC 
COORDINATES 

In this section, we see how the probability distribution 
function of the work W exerted on a system can be eval- 
uated by considering the joint probability distribution of 
W and the microscopic state of the system. This equa- 
tion was derived in refs. |7[ ^ (see also [|| ) . Let us first 
consider a system whose microscopic state i can take on 
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a finite number of values. To each such state is assigned 
an energy value Hi (fj.) , where /i is a parameter which is 
manipulated according to some protocol fi{t), starting at 
t = 0. We assume that the evolution of the system is 
described by a markovian stochastic process: given, for 
all pairs the transition rate fcij(t) from state j to 

state i at time t, the system satisfies the set of differential 
equations 

^=J2 [hAt)Pj{t)~k,MMt)], (2) 

where Pi{t) is the probability that the system is found 
at state i at time t. Let p1'^{^) represent the equilibrium 
distribution corresponding to a given value of /i. We have 

-/3//i(A.) 

PTit^) = (3) 

where — — e^^^i^ is the partition func- 

tion corresponding to the value fi of the parameter, and 
the corresponding free energy. We require that the 
transition rates kij{t) are compatible with the equilib- 



rium distribution p^'^{fi), i.e., that, for any i, 

J2 - fc,.(OPr(M(i))] - 0. (4) 

We assume that the system is at equilibrium at t — 0, 
and therefore, that Pi{t) satisfies the initial condition 

p,{t^O)^j^^if,iO)). (5) 

As pointed out in ref. |7|, the function Pi{t) does not pro- 
vide sufficient information on the work performed on the 
system during the manipulation process. We can how- 
ever consider the joint probability distribution ^i{W,t) 
that the system is found in state i, having received a work 
W, at time t. If the system is in the state i at time t, the 
infinitesimal work dWi done on it in the interval dt reads 

We have thus 



(7) 



r 



The last equality is obtained by substituting the expres- 
sion for 5Wi given in eq. (^) and by taking the first order 
expansion in 5t of the rhs. We are now able to write the 
set of differential equations which describe the distribu- 
tion functions ^i{W^t) 

dt 



(8) 



The joint probability distribution satisfies the 

initial condition 



<^,{w,Q)^5{w)pr'{m) 



(9) 



We are interested in the state-independent work prob- 
ability distribution P{W, t) defined by 



(10) 



It is convenient to introduce the generating function of 
with respect to the work distribution, defined by 



(11) 



(Notice that we adopt here, for later convenience, the 
opposite sign convention with respect to that adopted in 
ref. H].) We assume that ^i(W^t) vanishes fast enough, 
as \W\ — + oo, for ^'i(A, t) to exist for any A. The function 
'^i satisfies the initial condition 



exp[-/3ff,(^(0))] 



Z 



A'(O) 



(12) 



and evolves according to the differential equation 
9t*,(A,i)= j dW e^"^ dt^,{W,t) 



I 



= / dVP' e- 



9/1 



d_d^ 



^,{X,t). (13) 



Exploiting eq. (Q), it is easy to verify that if A = — /3, 
for any i at any time t, the solution of eq. (ttSf), with the 
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initial condition (|T^), reads 

_ _ cq 



>(0) 



V(o) 



pr(MW)- (14) 



We can thus straightforwardly verify the Jarzynski equal- 
ity: 



m(o) 



It is thus possible, in principle, to evaluate the probability 
distribution function of the work W by solving the equa- 
tions (||) or (|l3|) for all the microscopic states i. This 
approach has been implemented in ref. Q for a simple 
model of a biopolymer. 



II. COLLECTIVE VARIABLES 

The approach discussed in the previous section be- 
comes quickly unwieldy as the complexity of the system 
increases: the dimension of the system is equal to 
the number of microscopic states of the system. Clearly 
the system phase space must be sufficiently small for this 
protocol to be carried out, as in the case discussed in 0. 
In all the other cases, where the system considered is 
characterized by a large number of degrees of freedom, 
one usually introduces some collective variables, and an 
effective free energy, in order to reduce the complexity of 
the problem. The assumption underlying this approach 
is that the system reaches on a comparatively short time 
scale a quasiequilibrium state constrained by the instan- 
taneous value of the collective coordinate. Thus, on the 
the time scale of the experiment, the state of the sys- 
tem can be well summarized by the collective coordinate, 
with the corresponding free energy playing the role of the 
hamiltonian. 

Thus, we consider in the following a system character- 
ized by a generic equilibrium free energy function JF^(M), 
where /i is again the parameter which is manipulated, and 
M is some collective (mean-field) variable. (We shall con- 
sider in the following the case in which M is a scalar, but 
the analysis holds also if M is a collection of real vari- 
ables.) We assume that the system dynamics is stochastic 
and markovian: let P{M, t) denote the probability distri- 
bution function of the variable M at time t, then its time 
evolution will be described by the differential equation 



dP 

'dt 



HP. 



(16) 



where H is a. differential operator which depends on the 
parameter /i. We require that the operator Ti. is com- 
patible with the equilibrium distribution function of the 
system, i.e., that the relation 



He 







(17) 



holds for any value of 

The developments which follow were first obtained in 
ref. 1^ for a collection of noninteracting spins. 

We will consider a general mean-field system, described 
by a collective variable M and a generic free energy func- 
tion !F^{M). (The derivation can be easily generalized to 
the case in which M has more than one component.) The 
work done on a system during the manipulation, along a 
given stochastic trajectory M{t), is given by 



W ■ 



At' ix{t') 



dT^{M{t')) 
djj, 



(18) 



Using the same arguments as for the discrete case, one 
finds that the time evolution of the joint probability dis- 
tribution ^{M,W,t) of M and W is described by the 
differential equation 



(19) 



It can be easily shown that the solution of eq. (|T^) satis- 
fies the Jarzynski equality (0) identically |H . 

Equation ( [l9| ) becomes much easier to treat if one in- 
troduces the generating function ^'(Af, A,i) for the work 
distribution: 



Equation (O) becomes thus 



at ofi 



with the initial condition 



*(M,A,0) 



-/3^^(0)(M) 
Z^,(0) 



(20) 



(21) 



(22) 



These equations are exact for a collection of free spins, 
or for a mean-field Ising model. The partial differential 
equation (21) replaces the 2^ ordinary differential equa- 
tions (|l|), with i e {-1,-1-1}^, that one would obtain 
without the use of the collective coordinate M. 

We now derive a path integral representation of the 
solution of eq. (pl|), taking for the differential operator 
H the expression 



H 



E 

fc=0 



{9kiM)-}. 



(23) 
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(The coefficients gk{M) also depend on /z, but tliis de- 
pendence is understood to ligfiten tfie notation.) Let us 
introduce tfie generating function of ^'(M, X,t): 



(24) 



In tfie continuum limit, eq. ( |30| ) becomes 

/„M(tf)=M 
dA/o / VjVM 
Jm{0)=Mo 

exp{5[7,M]} *(Afo,A,0), (32) 



where 



Multiplying both sides of eq. ( pi| ) by exp(— 7Af), and 
integrating over M, we obtain 



5[7, A/] = [ dt C{t). 
Jo 



(33) 



OMj, = j dM e-'^^' (n + A/i d^T^^ ^ The "lagrangian" C is given by 



dM e"^*^ 
dAT e"'^*^ 



Qk 



*. (25) 



Then the fimction $1(7, A, i) satisfies 

f)(7, X,t + 6t) = J dM e-^^^ 

{l + St[H{-i,M) + \iid^F^,]}^, (26) 
where the function 7i(7, Af) is defined as 

H(7,M)=^7'=5fc(A//). (27) 

Given ri(7, A, i), we can evaluate ^(Af, A, t) from the ex- 
pression 

/+ioc 1 
^e^*^l)(7,A,i). (28) 

(In the following, we shall understand the integration lim- 
its on 7.) We obtain therefore 

*(A/, \t + 6t)= j dM' eT(^^-^^') 

{l + 5t [H(7, M') + Xfi *(Af', A, t) 

— [ dM' e^^^^"*^')+''*['"*'^'*'^'^+^''"^"-^"] 
27ri7 

X ^'(Ar', A,t). (29) 



Iterating, we obtain 



^{M,\t + Nt5t)= /dMo /"jj- 



d7idA/i 



27ri 



xexp{5[7,A/]} *(Afo,A,0), 
where the "action" 5 [7, M] is given by 



5{M - Mt) 
(30) 



Nt 



5[7,Af] = ^{7,(A/,-A£,_i) (31) 
1=1 

+ [H(7., M.) + \^i^^T^^t^){M,)\ ] . 



'zr 

C{t) = (7A/ + H(7,A/)-f A/i^ 



7(i),Af(t),,<(t) 



(34) 

Let N indicate the size of the system, and let us define 
the "intensive quantity" m = M/N . We can thus define, 
in the thermodynamic limit N go, m = const., the 
densities 



H{l,m) 



= lim 

N^oo 

= lim 



N ' 
n{-f,Nm) 
N 



(35) 
(36) 



The Lagrangian density "per spin" then reads 

£{t)^ lim ^ ^jm + H{j,m) + Xfi^. (37) 

Af— >oo iV Ufl 

In this way, the path integral appearing in eq. ( ^ ) as- 
sumes a form suitable for a saddle-point approximation 
for large system sizes N, as pointed out in [|[ ||. The 
parameter N plays a role akin to the inverse of Planck's 
constant h in the quasiclassical approximation of Feyn- 
man's path integral for quantum amplitudes jl^ . The 
result is the leading term in an asymptotic expansion in 
powers of N~^, which corresponds to the mean-field solu- 
tion of a statistical model. In ref. Q it was shown that the 
approximation works well for free spins. In ref. it was 
shown that for a mean-field spin system above the phase 
transition the approximation works rather well for sys- 
tem sizes N of the order of 10 and larger, but deteriorates 
as the transition is approached. It would be interesting 
to investigate in full the behavior of a finite-size system, 
in a situation when the corresponding infinite-sized sys- 
tem exhibits a phase transition. For a sufficiently fast 
manipulation protocol, in a large but finite system, the 
probability that a fluctuation overcoming the free energy 
barrier spontaneously arises should be very small. We 
expect therefore that the results of the infinite-size limit 
should hold better for faster protocols than for slower 
ones. These issues will be dealt with in future work. 

In the leading approximation, the path integral in 
eq. ( |3^ ) is dominated by the classical path (7c (0: '^c(O); 
solution of the differential equations 



SS 

SS 

Sm(t) 



= 



dH 

— + A/i 
dm dmdfj. 



(38) 



(39) 
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We shall now see that the requirement that the system is 
in equilibrium before the manipulation starts, imposes an 
initial condition on these equations. In order to evaluate 
the integral over Mq in eq. ( |3^ ) with the saddle-point 
method, we note that ^'(-M, A, 0) appearing on its rhs, 
is given by eq. (^). Furthermore, from the definition of 
i{t), eq. (||), it follows that 

dte{t) = TOtf7tf-mo7o+ / dt [-7m + H + A/ia^/^] 
)o Jo 

_ _ (40) 

Thus, substituting eq. ( [iol) into (p2|), and taking the 
derivative with respect to niQ — AIq/N, we obtain the 
saddle-point condition 



Thus, the equation of motions (38) and (39) have to be 
solved with the initial and the final conditions (41) and 
(|^: let (7*(0: m*{t)) denote the solution of eqs. (|38|j3|) 
satisfying these conditions. For each value of A, taking 
into account the initial value condition (|2^), the follow- 
ing saddle point estimation for r(A,tf) is obtained by 
eq. (M: 



T{X,tt) oc 



exp{-iV5(A)} 



where 



7(A)=/3/^„(m5) 



dt 4*(i). 



(48) 



(49) 



7(t=0) = -/3 



dm 



(41) 



t=o 



In this way one can devise a strategy to evaluate 
^'(M, A, if) for a given manipulation protocol fJ-(t), when 
the system size N is large enough. One has to solve 
the classical evolution equations (t 



In this equation, £*{t) is £{t) evaluated along the classical 
path (jrjt), m^(t)). In order to evaluate the integral on 
the rhs of eq. (|44|), we use the saddle point method again, 
and obtain 

P{Nw,tf) ^J\feyii>{-N[X*{w)w + g{X*{w))]}, (50) 



boundary condition: namely, eq. (O) should be imposed 
at t = 0, and the condition Nm{t{) = M should be im- 
posed at the final time tf. Once the relevant classical path 
(7c(t), 7T!,c(i)) has been evaluated, one can obtain the ac- 
tion density s[7c, rnd = limAr^oo Sljc, Nrnd/N from the 
expression 



M) with a two-point ^^ere A* (w) is the solution of 



.9'(A*) 



-w, 



(51) 



S be 'Tic 



dt l{t). 



(42) 



Then, taking into account the initial condition (22), 
we obtain the following asymptotic expression for 
'^{Nm,\ti): 

^'(A^m,A,tf) cxexp{Af [s[7e,me] -/3/^(o)(me(i=0))]}. 

(43) 

However, we are essentially interested in the state- 
independent work probability distribution 



and A/" is a normalization constant. Notice that the sad- 
dle point estimate for P{W, t{) obtained in this way, 
implies that the distribution becomes more and more 
sharply peaked around its maximum value as A'^ — > cx). 
This is compatible with the expectation that the work 
fluctuations becomes relatively smaller as the size of the 
system increases, and in the limit N — > 00, which can be 
thought as the limit of a macroscopic system, no work 
fluctuations are observed, and the work done on the sys- 
tem during the manipulation takes one single value, cor- 
responding to the most probable value of P{W,ti). 

In ref. [gj we showed that the JE is identically satisfied 
at the level of classical paths. For com pleteness, this 
derivation is reproduced in the Appendix . 



r(A,tf), 



where we have defined 



r{\,tf) = J dM *(Af, A,if). 



(44) 



(45) 



III. 



A MEAN-FIELD SYSTEM WITH 
LANGEVIN DYNAMICS 



We shall now see that evaluating r(A, ti) identifies a well- 
defined boundary condition on 7c (if). We have indeed 



r(A,tf) 



M{tt)=M 

dMdMo I V-fVM 

M{0)=Mo 



X exp 



N / dt£{t) 



*(Afo,A,0). 



(46) 



We wish to discuss a few properties of the work dis- 
tribution obtained by the present method by consider- 
ing a definite example. The case of free Ising spins has 
been considered (within a slightly different formalism) in 
ref. 1^. We shall return to it in sec. We thus take an 
Ising-like system with mean-field interaction, with free 
energy 

T{M) = -^M^ - hM - TS{M), (52) 

where where S{M) is the usual entropy for an Ising para- 
magnet. 



In order to evaluate the integral over M with the saddle 
point method, we notice that, upon derivation of the rhs 
of eq. (^) with respect to rriti, '^'^ obtain the condition 



S{M) = -fcE 



7f = 7(*f ) = 0- 



(47) 



N + M\^ fN + M 
log 



log 



(53) 
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expressed as a function of the continuous variable M. We 
assume that the system evolves according to Langevin 
dynamics. The corresponding Fokker-Planck differential 
operator reads 



n- =uJoN 



d 
'dM 



dM 



leading to the hamiltonian 



dm 



(54) 



(55) 



where the free energy density f{m) is given by 



/(m) = ——w? — hm + l3 ^ 



1 + m\ f 1 + in 
log 



1 — ?7i\, f \ — m 
log 



(56) 



The stochastic process described by this operator can 
be simulated by integrating the corresponding Langevin 
equation, using the Heun algorithm : for each realiza- 
tion of the process, the work W done on the system can 
be evaluated. The resulting histogram of w represents 
an estimate of the work probability distribution. This 
estimated distribution can be then compared with the 
expected distribution (valid asymptotically for N — + oo) 
obtained by the classical paths. 

We consider the case where the system is subject to 
the external manipulation of the magnetic field h{t)^ ac- 
cording to the simple protocol 



h{t)^hQ + {hi-h^)- 



0<t<t{. 



The equations of motion ( 

dH 

o-f 



391) become 



dl 
dm 

7 - — + Xfi-^^ 
dm dmdfi 



(57) 



(58) 



dm"^ 



7 - \h, (59) 



In the following we will take (3 = 1. In figure ^ we 
consider the case where the system is above the critical 
temperature, i.e., [3 J < 1. In this case, as expected, the 
peak of the distribution moves towards the value of the 
work done on the system along a reversible trajectory 
Wiov = (^Z = 0, as the transformation becomes slower. 
But the most important indication emerging from such 
a figure, is that the JE cannot be applied to obtain an 
independent estimation of the free energy difference be- 
tween the final and initial states of the transformation, 
if N is too large. In fact, we plot in the same figure, the 
quantity P{w) defined as 



P{w) = exp [~PNw] P{w), 



(60) 



on the one hand we find / dwP{w) = exp [— /3AF] = 1 
as predicted by the JE, while on the other hand the his- 
togram obtained by the simulations exhibits no point (no 
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FIG. 1: Results for the system described bv the differential 
operator (|^ with equihbrium free energy (^ij), manipulated 
according to the protocol (^7|), with J = 0.5, ho = ~hi = 
— 1, and {a)t{ = 2, {b)tf = 4 . Continuous line: probability 
density P{w) of the work "per spin" w = W/N, with A'^ = 100. 
The histogram of the work is obtained by 10000 simulations of 
the process, see text. Dotted line: P(w) as given by eq. (|60[), 
whose integral verifies the Jarzynski equality. Vertical line: 
Thermodynamic value of the work Wrcv = AF/N. 



realization of the process) with w < = Wrov Thus the 
work distribution obtained by the simulation of the pro- 
cess cannot reliably be used for estimating AF. This is a 
typical example of how the lack of knowledge of the tails 
of the work distributions in micro-manipulations experi- 
ments hinders the possibility of using eq. (|^) to evaluate 
free energy differences. 

We now consider a system below the transition tem- 
perature, i.e., for pj > 1. In figure ^ the work probabil- 
ity distribution obtained by the theory here discussed, is 
plotted for J = 1.1, ho = ~hi = —1, and for two values 
of the final time t{. In the same figure, the probabil- 
ity distribution obtained by simulations is also plotted. 
As for the case /3J < 1 (fig. |l|), the JE is satisfied, i.e., 
(exp {—PW)) = 1, there is a good agreement between the 
theory and the histograms obtained by simulations. But 
also in this case, such simulations cannot be used for esti- 
mating AF, since the histograms exhibits no point with 

W < Wrev 
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Since the amplitude of work fluctuations is expected 
to be relatively large in small system, we calculate now 
the work probability distribution for smaller systems and 
compare them with the results of simulations. First, we 
consider the case N = 10, fig]|: it can be seen that 
the the histogram of the work obtained by simulations is 
closer to the thermodynamic value of the work Wrcv = 0, 
than the distribution function obtained by the theory dis- 
cussed in the present paper. Indeed, since P{Nw,t), as 
given by eq. (pG), is exact only in the limit TV — > oo, that 
expression fails to describe the actual work distribution 
for small N. Furthermore, even for N — 10, there are 
few points in the histogram with w < Wrev, and thus no 
reliable estimate of AF can be obtained from the simu- 
lations. 

We further decrease the value of N and take N — 
2, see fig. In this case the agreement of the his- 
togram with the theoretical curve is worse than the 




FIG. 2: Results for the system described by the differen- 
tial operator (jH^) with equilibrium free energy (^^, manip- 
ulated according to the protocol (^7[), with Jo = Ji = 1.1, 
ho = —hi = —1, and (a)tf = 2, (h)ti = 4. Continuous line: 
probability density P{w) of the work "per spin" w — W/N, 
with A'' = 100. The histogram of the work is obtained by 
10000 simulations of the process, see text. Dotted line: P{'w) 
as given by eq. (^o|), whose integral verifies the Jarzynski 
equality. Vertical line: Thermodynamic value of the work 

Wrev = AF/N. 
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W 

FIG. 3: Results for the system described by the differential 
operator with equilibrium free energy (^2|), manipulated 
according to the protocol (^7|), with Jo = Ji = 0.5, ho = 
—hi = —1, and tf — 2. Continuous line: probability density 
P{w) of the work "per spin" w = W/N, with iV = 10. The 
histogram of the work is obtained by 10000 simulations of the 
process, see text. Dotted line: P(w) as given by eq. (|6o|), 
whose integral verifies the Jarzynski equality. Vertical line: 
Thermodynamic value of the work Wrov ~ AF/N. 

case TV — 10, as expected. But the small size of the 
system entails a broader work distribution, and thus 
enables a sufficient sampling of trajectories with w < 
w^ev In the same figure, the histogram of the distribu- 
tion exp [-/SNw] P{w) is plotted: from this histogram 
we obtain the estimate for the free energy difference 
A/exp - -N-^kBT\n[exp{-pNw) P{Nw)l^p, where 
[■ ■ -Icxp mean over all realizations of the process. 

We obtain A/exp — 0.015, against a theoretical value of 
Wrov = A/ = 0, and a most probable value of the work 
WniD - 0.6. 



IV. PATH THERMODYNAMICS 

The work distribution can be interpreted in terms of 
path thermodynamics, as first suggested in ref. [^. In- 
deed, g{X) = — hmAT^oo logr(A, tf)/N plays the role of a 
path Gibbs free energy. Thus 

</.H = - lim ^\ogP{Nw,tf), (61) 

plays the role of the corresponding Helmholtz free energy. 
The two functions are related by a Legendre transforma- 
tion: 

(t){w) = inf (.g(A) -|- Xw) 

^ g{X*{w))+X*{w)w, (62) 

where X*{w) is the solution of eq. (^l|). Thus A and w ap- 
pear like thermodynamically conjugate variables. Notice 
that if {X,w*{X)) are a pair of mutually conjugate vari- 
ables, then w*(A) is a monotonically increasing function 
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we obtain straightforwardly 

(j){m, w) ~ inf (cj(7, A) + 7m + Xw) , (65) 

where 01(7, A) is defined in terms of Q{X,j,t), which we 
have defined in eq. (24), by 



u.(7,A) = - lim -logr!(7,A,if). (66) 

One may notice that the 7 appearing in this equation 
may be identified with 7f = 7(if). 



FIG. 4: Results for the system described by the differen- 
tial operator (jH^) with equilibrium free energy (^^, manip- 
ulated according to the protocol (^7[), with Jo = Ji = 0.5, 
ho — —hi — —1, and tf = 2. Continuous line: probabil- 
ity density P{'w) of the work "per spin" w = W/N, with 
N — 2. The histogram of the work is obtained by 10000 
simulations of the process, see text, and is plotted with 
full line. Dotted line with black diamonds: histogram of 
P{'w) = exp{—PNw)P{w). Vertical line: Thermodynamic 
value of the work Wicv ~ AF/N. 



V. LARGE FLUCTUATIONS AND 
EXPONENTIAL TAILS 

It was suggested in ref. ||^, on the basis of numerical 
evidence, that, for slow protocols, the work distribution 
exhibits exponential tails. Here we discuss this intrigu- 
ing question. From eq. (|6T| ) we see that if P{Nw,ti) oc 
exp(— A^Aoit;) in some interval W- < w < W-^-, one has 



0(w) — Xqw + const. 



(67) 



of A. Indeed the relation between 0(w) and A reads 

(/)'(w*(A)) = A. (63) 

It is clear that the most probable value of the work w^p 
corresponds to the value A = 0. 

In ref. is taken to play the role of the internal 

energy, and thus —(j){'w) that of the entropy. Therefore 
A = (/''(w) can be considered as an inverse temperature. 
We have preferred to draw the analogy with more familiar 
functions. 

Indeed, one can generalize this point of view by go- 
ing back to the joint probability distribution function 
<^{M,W,t). If we define 



(j){m,'w) — — lim log^{Nm, NWjtf), 

N^oo 



(64) 



in the same interval. A linear behavior in the Helmholtz 
free energy is the signature of a first-order phase transi- 
tion. In the corresponding Gibbs free energy one has an 
angular point, i.e., a point Aq in which 



lim 



9'{X)=w.; l\mg'{X}^w+. (68) 

A^A+ 



Thus a horizontal plateau in a plot of A* vs. w corre- 
sponds to an exponential tail in P{Nw,tf). 

We shall now follow ref. p , by considering a system of 
N non interacting spins ct; = ±1, evolving according to 
the Glauber dynamics. The collective coordinate M is 
the total magnetization M = '^ii discrete variable) 
and the role of /i is played by the magnetic field h. The 
system evolves according to the master equation 



dP 
'dt 



N + M + 2 



P(M + 2,<) 



N + M 



PiM, t) 



N - M + 2 



N - M 



P{M,t) 



(69) 



where the spin flip rates p^'^ are given by 



(70) 



P{M -2,t)- 

I 

spin flip. In this case, the free energy Th{M) reads 

ThiM) = -hM ~TS{M), (71) 



where S{M) is given by eq. ( p4[ ) as a function of M. 
in which ujo{h) is a microscopic "attempt frequency" for We can make the connection with our formalism by 
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momentarily considering M as a continuous variable, and The master equation (|69|) then assumes the form 
by expressing the shift operator 



T±/(Af) = /(Af±2), 
in the following way 

T±-e±27&. 



(72) 



dP 
'dt 



= HP, 



(73) where the differential operator Ti. is given by 
I 



n- = 



N + M 



N 



(74) 



r 



where it is understood that the derivative on AI acts 
on all instances of M it finds on its right. Then the 
hamiltonian as given by eqs. ( p7| ) and (|3^) has the 
form 



H 



(e^^ - 1) 



1 + m 



-P 



— 1) z p 



(75) 



Equations ( p^ 
classical path: 



2 " ' '2 
yield the equations of motion for the 



"%T(1 - m) - e2V(l + m), 



(76) 



7=2 [("^'^ - 1) - (e"'^ - 1) p1 - (77) 

A different and more complicated approach, used in 
ref. 1^, leads to the same results. 

We shall suppose that the applied magnetic field h is 
manipulated according to the simple protocol (pTj) 



h{t) ^ho + (hi - ho) 



0<t<tf. 



We also suppose that uJo{h) = LUo/{e^^ + e '''') so that 
the functions , , are explicitly given by 



pHt) = UJQ 

p^{t) 



UJo 



Ql3h{t) Q-I3h{t) 
,J3h(t) Q~ph{t) 



(78) 

(79) 



where ujq is a constant. 

Let us now consider the quasi-static limit /i — > 0, with 
Xh ^ K = const. It is then possible to neglect the Ihs of 
eqs. (^60), yielding 

m = tanh (/3/i - 27) , (80) 
2k = {e^'^ - 1) - (e"^'^ - 1) . (81) 

Combining these equations, one obtains an expression for 
m as a function of h: 



sinh(/3ft,) — 2KC0sh(/3/i) 
^Jl + [smh{(3h) - 2KCOsh(/3/i)]^ 



(82) 



Thus rric depends on t via h, in terms of this equation. 
It also depends on the parameter k. One can check that 
mc{t, k) exhibits an extremum as a function of t in the 
interval [0,if], if |k| > Kc = 1/2, otherwise it is strictly 
monotonic. In order to discuss an explicit example, we 
set p = 1, hi = —ho = 10. In fig. ||, the function mdt, k), 
as given by eq. (|2|) is plotted for three different values of 
the parameter k: the function clearly exhibits a different 
behavior for k < Kr and n > Kr- Thus, as k becomes 




-0.5 - 



200 



FIG. 5: Plot of mc{t, k) as a function of t, as given by eq. (p2[), 
for three values of the parameter k. The external field is ma- 
nipulated according to eq. (^^, with t{ = 100. The function 
is monotonic for k < Kc. 

greater than its critical value Kc, we expect a singular 
behavior of the curve (k, w(k)), where 



w{k,) 



dt h{t)mc{t, k), 



(83) 



and mc(t, k) is given by eq. wAL Evaluating w{k) we 
obtain the curve plotted in fig. |6f one can see that it 
exhibits, for | pronounced minimum in An /Aw 

rather than a horizontal plateau. 

The simplicity of the system allows us to check this 
prediction by directly solving the equations ([T^ ) for the 
generating function "^^{X^t), a = ±1, for the transition 
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FIG. 6; Plot of K as a function of w in the quasistatic limit, 
as given by eq. (fea). The line is a guide to the eye. 



rates given by eqs. (|78|j79| ) . One thus obtains the function 
ri(A,it) for the single spin, from the expression 



(84) 



=±1 



Since g{X) ~ — log [FifA, tf)], we can obtain the curve 
{w,X*{w)) from eq. (pl[), and compare it with the 
predicted curve for the quasi-static limit, as given by 
eq. (HI). Such a comparison is shown in fig. 0: as the 
value of /i = {hi — ho)/tf decreases the agreement be- 
tween the theory and the curve predicted by eq. (^) 
improves. 
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1 
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FIG. 7; Plot of At = A/i as a function of w, for different values 
oi h — {hi — hO)/t{. Thick full line: k{'w) in the quasi-static 
limit, as obtained by eq. (|83|). The other curves are obtained 
from the evolution equations (|l|). Dashed line tf = 20 (/i = 
1), dotted line tc = 200 {h = 0.1), the curve with tc = 2000 
{h = 0.01) is practically indistinguishable from the quasi- 
static limit. 

We checked that this behavior depends on the details of 
the dynamics by considering the same paramagnetic sys- 
tem, but evolving by a Langevin rather than a Glauber 



dynamics, with the method reported in section II] 
shown in fig. 



As 



the corresponding (w, k) curve exhibits 
no plateau, and therefore there are no exponential tails 
in the work distribution. These results are confirmed by 
a detailed analysis of the quasi-static limit. 



3 r- 
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1 

!^ - 
-1 - 





w 



FIG. 8: Plot of K as a function of w for the paramagnet 
evolving according to a Langevin dynamics. Manipulation 
protocol (p7[), with ho — —hi — —5, tf = 1000. 

Let us now turn again to the ferromagnetic mean-field 
system with Langevin dynamics. 

In the top panel of figure ^, we plot m* as a function of 
t for different values of A, obtained by numerical solution 
of eqs. ( |5^ , |59| ), for J = 0.5 and tf = 2. It can be seen 
that the shape of m*{t) varies continuously as A is varied. 
Accordingly there is no horizontal plateau in the A* vs. 
w plot, implicitly defined by eq. (^l|), as shown in the 
bottom panel of fig. ^. The same behavior is obtained by 
varying U and implementing a slower or a faster protocol 
(data not shown). According to the above discussion, 
the work distribution exhibits no exponential tails for 
the case (3 J < 1. 

We now investigate whether a different behavior can 
appear when the system is manipulated across the 
symmetry-breaking transition it exhibits at /3J = 1. Let 
us look at the behavior of the classical path m* (t. A) cor- 
responding to the 7f = boundary condition, both above 
{(3 J < 1) and below (/3J > 1) the transition. 

In the top panel of figure |lO|, we plot m* as a function of 
t for different values of A, obtained by numerical solution 
of eqs. (^, ^), for J = 1.1 and tf = 2: we observe no 
discontinuity in m*(i. A) as A is varied, and thus w is a 
continuous function of A, see figure |l^ bottom panel. 

We now consider a faster protocol, tf = 0.2, with the 
same value of J and Hq: the results are plotted in figure 
|l|. One can clearly see that TO*(t, A) exhibits a discon- 
tinuity for A = 0.5, jumping from negative to positive 
values. Accordingly, w(A*) exhibits a discontinuity at 
A* 
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0.5, as shown in the bottom panel of fig. 
If we now evaluate the path Helmholtz free energy 4>{w) 
from eq. ( |6^ ) , we obtain the results shown in fig. |l^. As 
discussed in section IV, 4'{w) should be obtained by a lin- 



11 




-1 ' — • 1 

0.5 1 1.5 2 

t 




-2 -1.5 -1 -0.5 0.5 1 1.5 2 
W 



FIG. 9; Top: plot of m* as a function of t for different values 
of A, with J = 0.5, ho = —hi = —1, and ti — 2. The values of 
A vary between A = — 5 (bottom curve) and A = 5 (top curve), 
with a step AA — 0.2. Bottom: plot of A* as a function of w, 
as defined by eq. (fell). 
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FIG. 10: Top: plot of m* as a function of t for different values 
of A, with J — 1.1, ho = —hi — —1, and t{ — 2. The values of 
A vary between A — —5 (bottom curve) and A = 5 (top curve), 
with a step AA = 0.2. Bottom: plot of A* as a function of w, 
as defined by eq. (fell). 



ear interpolation between [w-i-, (j){w^)) and (w-, 0(ii'_)), 
where 'w± are the values of w either side of the discon- 
tinuity. This corresponds to an exponential tail in the 
distribution of the work. In this case, the existence of an 
equiUbrium phase transition shows up as a path phase 
transition, i.e., an exponential tail, provided that the ma- 
nipulation protocol is fast enough. The same behavior is 
obtained for ho = —hi = —0.1, tf = 0.2 (no discontinu- 
ity) and tf = 0.02 (discontinuity, data not shown). This 
last result suggests thus that the presence of exponential 
tails in the work probability distribution is due to a path 
"phase separation" , which is induced by a sufhciently fast 
manipulation protocol: inspection of fig. |ll| indicates that 
the trajectories ?n*(t. A) form two groups, as A is varied, 
and none of the trajectories belonging to each of the two 
groups crosses the line m — 0, differently from what hap- 
pens for a slower protocol, see fig. |l^. We checked that 
the resulting distribution (f>{w) satisfies the following re- 
lation, which is a consequence of Crooks' identity ||] and 
of the symmetry h{tf — t) = —h{t) satisfied by our pro- 
tocol: 

(l){w) - 4>{-w) = -f3w. (85) 



It would be interesting to see if such a "path phase tran- 
sition" takes place in more realistic models. 

VI. DISCUSSION 

In this work, we have examined the distribution of the 
work W exerted on a system which is manipulated out of 
equilibrium. We have first obtained its expression by con- 
sidering the joint distribution of the microscopic state of 
the system and of the work. The expression one obtains is 
in principle exact, but is amenable to numerical solution 
only for very simple systems. We have then considered a 
system whose quasiequilibrium state can be described by 
one (or more) collective variables, to which an effective 
free energy function is associated. The resulting equa- 
tion for the joint distribution of the collective variables 
and work is a partial differential equation which can in 
principle be numerically solved. However, we found that 
it is possible to explore a different direction. Indeed, fol- 
lowing ref. 1^ , one sees that one can express the solution 
to this equation as a path integral. In the limit of system 
size N going to infinity, the path integral is dominated by 
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FIG. 11: Top: plot of m* as a function of t for different values 
of A, with J = 1.1, ho — —hi — —1, and tf — 0.2. The values 
of A vary between A = —5 (bottom curve) and A = 5 (top 
curve), with a step AA = 0.2. Thick line: A = 0.5. Bottom: 
plot of A* as a function w, as defined by eq. (|5l|). 

the classical paths, which satisfy a "canonical" system of 
ordinary differential equations, with suitable boundary 
conditions. Building on this information, it is possible to 
estimate the work probability distribution function for 
large system size, in the form 

P{W) oc exp[-A^(/)(w)] , 

where w = W/N, and (/>(u') plays the role of a work free 
energy density, or of a function of large deviations. This 
quantity is obtained as a Legendre transform of g(X) as 
given by eq. (^. It is natural to interpret the relations 
between these quantities as corresponding to those be- 
tween the Helmholtz and the Gibbs free energy densities 
in thermodynamics. Within this picture, the parameter 
A can be viewed as the intensive field conjugated with 
the extensive variable w, which acts as an order param- 
eter for the single path. Thus, horizontal plateau in the 
A* vs. w plot indicates a first-order phase transition in 
the paths. In this case the work distribution exhibits an 
exponential tail in a given range of w, depending on the 
manipulation details. Our results suggest that the system 
exhibits such path "phase separation" for sufficiently fast 
manipulation protocols, and below the mean-field equi- 
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FIG. 12: Path Helmholtz free energy for the system described 
by the differential operator (^^, manipulated according to the 
protocol (E^, with J = 1.1, ho = —hi = —1, and tt = 0.2. 



librium transition temperature, whereas above it one can 
find only a marked infiection point in the A* vs. w plot, 
but not a horizontal plateau. 

The results we obtain are interesting in their own right, 
since they exhibit a number of nontrivial properties of 
the classical paths. However, their usefulness for assess- 
ing the feasibility of the use of the Jarzynski equality for 
the reconstruction of the equilibrium free energy land- 
scape can be a priori doubted. Indeed, the P{W) one 
obtains in this way is only asymptotically valid for large 
N, and in this case, the probability of observing, in an ac- 
tual experiment, a sufficient number of large fiuctuations 
to evaluate the Jarzynski average (|l|) with some confi- 
dence, is extremely small. We found however that the 
estimated distribution is not too far from the actual dis- 
tribution for system sizes as small as 2, at least when the 
manipulation protocol is not too fast and does not cross 
an equilibrium phase transition line In this case the 
JE can be successfully applied to the work distribution 
obtained by simulations: the estimate of the free energy 
difference differs little from the expected value. 

It is reasonable to expect, for our mean-field like sys- 
tems, that the existence of a first-order transitions could 
cause some problems. Formally, in the limit iV — > oo and 
for a manipulation protocol with a finite speed, the sys- 
tem would remain close to the free-energy minimum it 
finds itself in until it reaches the spinodal line. In a finite 
system, if the protocol is slow enough, the system can 
cross the free energy barrier and reach the real minimum 
in a finite time. We found that the classical paths are 
able to interpolate between the minima for slow enough 
protocols, whereas they tend to split in different phases 
for fast ones. Thus this effect takes place even for mean- 
field systems. 

It is possible to extend this work to more realistic sys- 
tems, provided that the basic assumption of the existence 
of relevant collective variables holds. One should also 
consider what information can be gathered by exploiting 
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other manipulation protocols. 



Let us now take the derivative of eq. (A.l) with respect 
to 7 at fixed fj,. We obtain 



Acknowledgments 

We thank F. Ritort for his interest in our work. This 
research was partially supported by MIUR-PRIN 2004. 



d'^ff_c dm 



(A.5) 



By multiplying both sides of eq. (A.4) by d f^/dm and 



substituting eq. (A.5), we obtain the following relation 



APPENDIX: DERIVATION OF THE JARZYNSKI 
EQUALITY FOR THE CLASSICAL PATHS 



dm? dj dm 



(A.6) 



We report here, for completeness, the derivation of the 
Jarzynski equality at the level of classical paths, obtained 
in ref. [||. We first show that, for A = — /?, the solution 
of the classical equations of motion (|3|,^) satisfy an 
equation analogous to (ETl) at all times, namely 



Q 



om 



(A.l) 



By multiplying both sides of eq. ( pT| ) by e ''^^ and inte- 
grating by parts over M one obtains 



dM H(7, M) e-P^.iM)-iM ^ 0^ 



(A.2) 



where ?i{'j,M) is given by (p6|). Evaluating this inte- 
gral by the saddle point method in the large N limit, we 
obtain 



iJ(7,m*)=0, 



(A.3) 



if 7 and m* are related by ( |A.1| ). By differentiating 
eq. (A.3) with respect to 7 at fixed fj, we obtain 



dH 

d^ 



dH 

dm 



dm* 
97 



0. 



(A.4) 



which holds when 7 and m are related by eq. (A.l). 
We can now evaluate the time derivative of the Ihs of 

Mj. We have 



eq. (A.l), when 7 and m satisfy eqs. 



-7-/3 



dm^ 



d'U ^ 
dmd^ 



(A.7) 



dH ^ dy, 

dm dmdjji 



dH 



(3 



d^. 



-11. 



dm^ 97 dmdfx ' 
The seco nd and the last term cancel out. Substituting 
eq. (A.6), we see that also the first and the third therm 
cancel out. Thus if 7 and m satisfy eqs. (p8|,p9|) at all 



times, and satisfy eq. (A.l) at a given time, they satisfy 
this last equation at any time. 

Thus, for A = — /3, the Lagrangian, evaluated along the 
classical path, is given by 



N 



7m ■ 



N 



-p^m 
om 



(A., 



where we have exploited eq. (A.l). Substituting this ex- 
pression in eq. (32) one recovers eq. (|l^ ) and the Jarzyn- 
ski equality. 
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